1. Data Preprocessing

1.1 Gloabl Confirmed Cases and Deaths Time-Series by Country

 ts.confirmed = covid19.data("ts-confirmed")
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
## --------------------------------------------------------------------------------
 ts.confirmed =  ts.confirmed %>%
   group_by(Country.Region) %>%
   select(-c(Province.State, Lat, Long)) %>%
   summarise_all(funs(sum)) %>%
   arrange(desc(.[[ncol(ts.confirmed)-3]]))
## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.
 head(ts.confirmed)
## # A tibble: 6 x 97
##   Country.Region `2020-01-22` `2020-01-23` `2020-01-24` `2020-01-25`
##   <fct>                 <int>        <int>        <int>        <int>
## 1 US                        1            1            2            2
## 2 Spain                     0            0            0            0
## 3 Italy                     0            0            0            0
## 4 France                    0            0            2            3
## 5 Germany                   0            0            0            0
## 6 United Kingdom            0            0            0            0
## # … with 92 more variables: `2020-01-26` <int>, `2020-01-27` <int>,
## #   `2020-01-28` <int>, `2020-01-29` <int>, `2020-01-30` <int>,
## #   `2020-01-31` <int>, `2020-02-01` <int>, `2020-02-02` <int>,
## #   `2020-02-03` <int>, `2020-02-04` <int>, `2020-02-05` <int>,
## #   `2020-02-06` <int>, `2020-02-07` <int>, `2020-02-08` <int>,
## #   `2020-02-09` <int>, `2020-02-10` <int>, `2020-02-11` <int>,
## #   `2020-02-12` <int>, `2020-02-13` <int>, `2020-02-14` <int>,
## #   `2020-02-15` <int>, `2020-02-16` <int>, `2020-02-17` <int>,
## #   `2020-02-18` <int>, `2020-02-19` <int>, `2020-02-20` <int>,
## #   `2020-02-21` <int>, `2020-02-22` <int>, `2020-02-23` <int>,
## #   `2020-02-24` <int>, `2020-02-25` <int>, `2020-02-26` <int>,
## #   `2020-02-27` <int>, `2020-02-28` <int>, `2020-02-29` <int>,
## #   `2020-03-01` <int>, `2020-03-02` <int>, `2020-03-03` <int>,
## #   `2020-03-04` <int>, `2020-03-05` <int>, `2020-03-06` <int>,
## #   `2020-03-07` <int>, `2020-03-08` <int>, `2020-03-09` <int>,
## #   `2020-03-10` <int>, `2020-03-11` <int>, `2020-03-12` <int>,
## #   `2020-03-13` <int>, `2020-03-14` <int>, `2020-03-15` <int>,
## #   `2020-03-16` <int>, `2020-03-17` <int>, `2020-03-18` <int>,
## #   `2020-03-19` <int>, `2020-03-20` <int>, `2020-03-21` <int>,
## #   `2020-03-22` <int>, `2020-03-23` <int>, `2020-03-24` <int>,
## #   `2020-03-25` <int>, `2020-03-26` <int>, `2020-03-27` <int>,
## #   `2020-03-28` <int>, `2020-03-29` <int>, `2020-03-30` <int>,
## #   `2020-03-31` <int>, `2020-04-01` <int>, `2020-04-02` <int>,
## #   `2020-04-03` <int>, `2020-04-04` <int>, `2020-04-05` <int>,
## #   `2020-04-06` <int>, `2020-04-07` <int>, `2020-04-08` <int>,
## #   `2020-04-09` <int>, `2020-04-10` <int>, `2020-04-11` <int>,
## #   `2020-04-12` <int>, `2020-04-13` <int>, `2020-04-14` <int>,
## #   `2020-04-15` <int>, `2020-04-16` <int>, `2020-04-17` <int>,
## #   `2020-04-18` <int>, `2020-04-19` <int>, `2020-04-20` <int>,
## #   `2020-04-21` <int>, `2020-04-22` <int>, `2020-04-23` <int>,
## #   `2020-04-24` <int>, `2020-04-25` <int>, `2020-04-26` <int>
ts.deaths = covid19.data("ts-deaths")
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
## --------------------------------------------------------------------------------
ts.deaths =  ts.deaths %>%
   group_by(Country.Region) %>%
   select(-c(Province.State, Lat, Long)) %>%
   summarise_all(funs(sum)) %>%
   arrange(desc(.[[ncol(ts.deaths)-3]]))
 head(ts.deaths)
## # A tibble: 6 x 97
##   Country.Region `2020-01-22` `2020-01-23` `2020-01-24` `2020-01-25`
##   <fct>                 <int>        <int>        <int>        <int>
## 1 US                        0            0            0            0
## 2 Italy                     0            0            0            0
## 3 Spain                     0            0            0            0
## 4 France                    0            0            0            0
## 5 United Kingdom            0            0            0            0
## 6 Belgium                   0            0            0            0
## # … with 92 more variables: `2020-01-26` <int>, `2020-01-27` <int>,
## #   `2020-01-28` <int>, `2020-01-29` <int>, `2020-01-30` <int>,
## #   `2020-01-31` <int>, `2020-02-01` <int>, `2020-02-02` <int>,
## #   `2020-02-03` <int>, `2020-02-04` <int>, `2020-02-05` <int>,
## #   `2020-02-06` <int>, `2020-02-07` <int>, `2020-02-08` <int>,
## #   `2020-02-09` <int>, `2020-02-10` <int>, `2020-02-11` <int>,
## #   `2020-02-12` <int>, `2020-02-13` <int>, `2020-02-14` <int>,
## #   `2020-02-15` <int>, `2020-02-16` <int>, `2020-02-17` <int>,
## #   `2020-02-18` <int>, `2020-02-19` <int>, `2020-02-20` <int>,
## #   `2020-02-21` <int>, `2020-02-22` <int>, `2020-02-23` <int>,
## #   `2020-02-24` <int>, `2020-02-25` <int>, `2020-02-26` <int>,
## #   `2020-02-27` <int>, `2020-02-28` <int>, `2020-02-29` <int>,
## #   `2020-03-01` <int>, `2020-03-02` <int>, `2020-03-03` <int>,
## #   `2020-03-04` <int>, `2020-03-05` <int>, `2020-03-06` <int>,
## #   `2020-03-07` <int>, `2020-03-08` <int>, `2020-03-09` <int>,
## #   `2020-03-10` <int>, `2020-03-11` <int>, `2020-03-12` <int>,
## #   `2020-03-13` <int>, `2020-03-14` <int>, `2020-03-15` <int>,
## #   `2020-03-16` <int>, `2020-03-17` <int>, `2020-03-18` <int>,
## #   `2020-03-19` <int>, `2020-03-20` <int>, `2020-03-21` <int>,
## #   `2020-03-22` <int>, `2020-03-23` <int>, `2020-03-24` <int>,
## #   `2020-03-25` <int>, `2020-03-26` <int>, `2020-03-27` <int>,
## #   `2020-03-28` <int>, `2020-03-29` <int>, `2020-03-30` <int>,
## #   `2020-03-31` <int>, `2020-04-01` <int>, `2020-04-02` <int>,
## #   `2020-04-03` <int>, `2020-04-04` <int>, `2020-04-05` <int>,
## #   `2020-04-06` <int>, `2020-04-07` <int>, `2020-04-08` <int>,
## #   `2020-04-09` <int>, `2020-04-10` <int>, `2020-04-11` <int>,
## #   `2020-04-12` <int>, `2020-04-13` <int>, `2020-04-14` <int>,
## #   `2020-04-15` <int>, `2020-04-16` <int>, `2020-04-17` <int>,
## #   `2020-04-18` <int>, `2020-04-19` <int>, `2020-04-20` <int>,
## #   `2020-04-21` <int>, `2020-04-22` <int>, `2020-04-23` <int>,
## #   `2020-04-24` <int>, `2020-04-25` <int>, `2020-04-26` <int>

1.2 U.S. Confirmed Cases and Deaths Time-Series By State

nCov = load_nCov2019()
df.nCov = nCov$province
ts.us = df.nCov %>% filter(country == "United States") %>% 
  select("Date" = time, "State" = province, "Confirmed" = cum_confirm, "Deaths" = cum_dead)
ts.us$State = sapply(ts.us$State, Caps)
head(ts.us)
##         Date            State Confirmed Deaths
## 1 2020-03-15 Washington State       769     42
## 2 2020-03-15   New York State       746      6
## 3 2020-03-15       California       431      6
## 4 2020-03-15    Massachusetts       164      0
## 5 2020-03-15          Florida       136      5
## 6 2020-03-15         Colorado       131      1

1.3 Global Aggregated Cases and Deaths by Country

agg.gloabl = covid19.data("aggregated")
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
agg.gloabl = agg.gloabl %>%
  select(Country_Region, Confirmed, Deaths) %>%
  group_by(Country_Region) %>%
  summarise_all(funs(sum)) %>%
  arrange(desc(Confirmed))
head(agg.gloabl)
## # A tibble: 6 x 3
##   Country_Region Confirmed Deaths
##   <fct>              <int>  <int>
## 1 US                938154  53755
## 2 Spain             223759  22902
## 3 Italy             195351  26384
## 4 France            161644  22648
## 5 Germany           156513   5877
## 6 United Kingdom    149569  20381

1.4 U.S. Aggregated Cases and Deaths By State

agg.us = covid19.data("aggregated")
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
agg.us = agg.us %>%
  filter(Country_Region == "US" & Province_State != "Recovered") %>%
  select(Province_State, Confirmed, Deaths) %>%
  group_by(Province_State) %>%
  summarise_all(funs(sum)) %>%
  arrange(desc(Confirmed))
head(agg.us)
## # A tibble: 6 x 3
##   Province_State Confirmed Deaths
##   <fct>              <int>  <int>
## 1 New York          282143  22009
## 2 New Jersey        105498   5914
## 3 Massachusetts      53348   2730
## 4 California         42368   1689
## 5 Illinois           41777   1875
## 6 Pennsylvania       41153   1793

1.5 U.S. Active Cases, Deaths and Recovered Time-Series

ts.us.all = covid19.data("ts-ALL") %>%
  filter(Country.Region == "US") %>%
  select(-c(Province.State, Country.Region, Long, Lat)) %>%
  group_by(status) %>%
  summarise_all(funs(sum))
## Data being read from JHU/CCSE repository
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Reading data from https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv
## Data retrieved on 2020-04-26 22:49:33 || Range of dates on data: 2020-01-22--2020-04-26 | Nbr of records: 264
## --------------------------------------------------------------------------------
## Data being read from JHU/CCSE repository
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Reading data from https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv
## Data retrieved on 2020-04-26 22:49:34 || Range of dates on data: 2020-01-22--2020-04-26 | Nbr of records: 264
## --------------------------------------------------------------------------------
## Data being read from JHU/CCSE repository
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Reading data from https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv
## Data retrieved on 2020-04-26 22:49:34 || Range of dates on data: 2020-01-22--2020-04-26 | Nbr of records: 250
## --------------------------------------------------------------------------------
ts.us.all = as_tibble(t(as.matrix(ts.us.all)))
## Warning: The `x` argument of `as_tibble.matrix()` must have column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
colnames(ts.us.all) = unlist(ts.us.all[1,])
ts.us.all = ts.us.all[-1,]
ts.us.all$confirmed = as.integer(ts.us.all$confirmed)
ts.us.all$death = as.integer(ts.us.all$death)
ts.us.all$recovered = as.integer(ts.us.all$recovered)
ts.us.all = ts.us.all %>% mutate(active = confirmed - death - recovered)
head(ts.us.all)
## # A tibble: 6 x 4
##   confirmed death recovered active
##       <int> <int>     <int>  <int>
## 1         1     0         0      1
## 2         1     0         0      1
## 3         2     0         0      2
## 4         2     0         0      2
## 5         5     0         0      5
## 6         5     0         0      5

2. Exploratory Analysis

2.1 World Map of Confirmed Cases and Deaths

world_map <- map_data("world")
agg.gloabl.map = agg.gloabl
agg.gloabl.map$Country_Region = as.character(agg.gloabl$Country_Region)
agg.gloabl.map$Country_Region[agg.gloabl.map$Country_Region == "US"] = "USA"
agg.gloabl.map$Country_Region[agg.gloabl.map$Country_Region == "United Kingdom"] = "UK"
agg.gloabl.map$Country_Region[agg.gloabl.map$Country_Region == "Korea, South"] = "South Korea"
agg.gloabl.map$Country_Region[agg.gloabl.map$Country_Region == "Taiwan*"] = "Taiwan"

agg.gloabl.map <- inner_join(agg.gloabl.map, world_map, by = c("Country_Region" = "region"))
ggplotly(
  ggplot(agg.gloabl.map, aes(long, lat, group = group, fill = Confirmed, text = paste0(Country_Region, "\n", "Confirmed Cases:", Confirmed))) +
    scale_fill_gradient(low = "gray99", high = "coral3") +
    geom_polygon(), tooltip = c("text")
         )
ggplotly(
  ggplot(agg.gloabl.map, aes(long, lat, group = group, fill = log(Confirmed), text = paste0(Country_Region, "\n", "Confirmed Cases:", Confirmed))) +
    scale_fill_gradient(low = "gray99", high = "coral3") +
    geom_polygon(), tooltip = c("text")
         )
ggplotly(
  ggplot(agg.gloabl.map, aes(long, lat, group = group, fill = Deaths, text = paste0(Country_Region, "\n", "Deaths:", Deaths))) +
    scale_fill_gradient(low = "gray99", high = "coral3") +
    geom_polygon(), tooltip = c("text")
         )
ggplotly(
  ggplot(agg.gloabl.map, aes(long, lat, group = group, fill = log(Deaths), text = paste0(Country_Region, "\n", "Deaths:", Deaths))) +
    scale_fill_gradient(low = "gray99", high = "coral3") +
    geom_polygon(), tooltip = c("text")
         )

2.2 US Map of Confirmed Cases and Deaths

states_map = map_data("state")
agg.us.map = agg.us %>% mutate(Province_State = tolower(Province_State))
agg.us.map = inner_join(agg.us.map, states_map, by = c("Province_State" = "region"))
ggplotly(
  ggplot(agg.us.map, aes(long, lat, group = group, fill = Confirmed, text = paste0(Province_State, "\n", "Confirmed Cases:", Confirmed))) +
    scale_fill_gradient(low = "gray99", high = "coral3") +
    geom_polygon(), tooltip = c("text")
         )
ggplotly(
  ggplot(agg.us.map, aes(long, lat, group = group, fill = log(Confirmed), text = paste0(Province_State, "\n", "Confirmed Cases:", Confirmed))) +
    scale_fill_gradient(low = "gray99", high = "coral3") +
    geom_polygon(), tooltip = c("text")
         )
ggplotly(
  ggplot(agg.us.map, aes(long, lat, group = group, fill = Deaths, text = paste0(sapply(Province_State, Caps), "\n", "Deaths:", Deaths))) +
    scale_fill_gradient(low = "gray99", high = "coral3") +
    geom_polygon(), tooltip = c("text")
         )
ggplotly(
  ggplot(agg.us.map, aes(long, lat, group = group, fill = log(Deaths), text = paste0(sapply(Province_State, Caps), "\n", "Deaths:", Deaths))) +
    scale_fill_gradient(low = "gray99", high = "coral3") +
    geom_polygon(), tooltip = c("text")
         )

2.3 Global Time-Series Plot (Top 10 Countries)

top10_countries = ts.confirmed$Country.Region[1:10]
ts.confirmed.top10 = ts.confirmed %>% filter(Country.Region %in% top10_countries)
df.global.ts = tibble()
n_days = 100
for(i in c(1:10)) {
  # get the variables
  Count = unlist(ts.confirmed.top10[i,-1])[unlist(ts.confirmed.top10[i,-1]) > n_days]
  Day = c(1:length(Count))
  Country = rep(top10_countries[i], length(Count))
  # append to the pre-specified table
  df.new = tibble(Country, Day, Count)
  df.global.ts = rbind(df.global.ts, df.new)
}
ggplotly(
  ggplot(df.global.ts, aes(Day, Count, group = Country, color = Country, 
                         text = paste0(Country, "\n", "Total Cases:", Count))) +
    geom_line() +
    ggtitle("Global Total Confirmed Cases (Top 10 Countries)") +
    xlab("Days After 100 Cases") +
    ylab("Total Confirmed Cases"), tooltip = c("text")
  )

2.4 U.S. Time-Series Plot (Top 10 States)

top10_states = ts.us %>% group_by(State) %>% summarise(max = max(Confirmed)) %>% 
  arrange(desc(max)) %>% select(State)
top10_states = top10_states$State[1:10]
ggplotly(
  ggplot(ts.us %>% filter(State %in% top10_states), aes(Date, Confirmed, group = State, color = State,
                                                        text = paste0(State, "\n", "Total Cases:", Confirmed))) +
    geom_line() +
    ggtitle("U.S. Total Confirmed Cases (Top 10 States)") +
    xlab("Date") +
    ylab("Total Confirmed Cases"), tooltip = c("text")
  )

2.5 World Trend

ts.world = colSums(ts.confirmed[,-1])
ts.world = tibble("Date" = names(ts.world), "Count" = ts.world)
ts.world$Date = as.Date(ts.world$Date)
ts.world$New = c(NA, diff(ts.world$Count))
coeff = tail(ts.world$Count, 1)/tail(ts.world$New, 1)
newColor = "black"
countColor = "darkorange"

moving_n = 10
ggplot(ts.world, aes(x=Date)) +
  geom_line(aes(y=New), color=newColor, alpha = 0.2) +
  geom_ma(aes(y=New), ma_fun = SMA, n = moving_n, color=newColor) +
  geom_bar(aes(y=Count/coeff), stat="identity", fill=countColor, color=countColor, alpha=.2) + 
  
  scale_y_continuous(
    
    # Features of the first axis
    name = "Number of New Cases",
    
    # Add a second axis and specify its features
    sec.axis = sec_axis(~.*coeff, name="Total Confirmed Cases")
  ) + 
  
  
  theme(
    axis.title.y = element_text(color = newColor, size=11),
    axis.title.y.right = element_text(color = countColor, size=11)
  )
## Warning: Removed 1 row(s) containing missing values (geom_path).

2.6 U.S. Trend

ts.us.confirmed = unlist(ts.confirmed[1, -1])
ts.us.confirmed = tibble("Date" = names(ts.us.confirmed), "Count" = ts.us.confirmed)
ts.us.confirmed$Date = as.Date(ts.us.confirmed$Date)
ts.us.confirmed$New = c(NA, diff(ts.us.confirmed$Count))
coeff = tail(ts.us.confirmed$Count, 1)/tail(ts.us.confirmed$New, 1)

ggplot(ts.us.confirmed, aes(x=Date)) +
  geom_line(aes(y=New), color=newColor, alpha = 0.2) +
  geom_ma(aes(y=New), ma_fun = SMA, n = moving_n, color=newColor) +
  geom_bar(aes(y=Count/coeff), stat="identity", fill=countColor, color=countColor, alpha=.2) + 
  
  scale_y_continuous(
    
    # Features of the first axis
    name = "Number of New Cases",
    
    # Add a second axis and specify its features
    sec.axis = sec_axis(~.*coeff, name="Total Confirmed Cases")
  ) + 
  
  
  theme(
    axis.title.y = element_text(color = newColor, size=11),
    axis.title.y.right = element_text(color = countColor, size=11)
  )
## Warning: Removed 1 row(s) containing missing values (geom_path).

3. U.S. COVID-19 Analysis and Prediction

3.1 U.S. Growth Rate

data = covid19.data("ts-confirmed")
## Data being read from JHU/CCSE repository
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Reading data from https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv
## Data retrieved on 2020-04-26 22:49:53 || Range of dates on data: 2020-01-22--2020-04-26 | Nbr of records: 264
## --------------------------------------------------------------------------------
growth.rate(data, geo.loc="US")
## Processing...  US
## Loading required package: pheatmap
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:PerformanceAnalytics':
## 
##     textplot
## The following object is masked from 'package:stats':
## 
##     lowess

## $Changes
##   geo.loc 2020-01-23 2020-01-24 2020-01-25 2020-01-26 2020-01-27 2020-01-28
## 1      US          0          1          0          3          0          0
##   2020-01-29 2020-01-30 2020-01-31 2020-02-01 2020-02-02 2020-02-03 2020-02-04
## 1          0          0          2          1          0          3          0
##   2020-02-05 2020-02-06 2020-02-07 2020-02-08 2020-02-09 2020-02-10 2020-02-11
## 1          0          0          0          0          0          0          1
##   2020-02-12 2020-02-13 2020-02-14 2020-02-15 2020-02-16 2020-02-17 2020-02-18
## 1          0          1          0          0          0          0          0
##   2020-02-19 2020-02-20 2020-02-21 2020-02-22 2020-02-23 2020-02-24 2020-02-25
## 1          0          0          2          0          0         36          0
##   2020-02-26 2020-02-27 2020-02-28 2020-02-29 2020-03-01 2020-03-02 2020-03-03
## 1          6          1          2          8          6         24         20
##   2020-03-04 2020-03-05 2020-03-06 2020-03-07 2020-03-08 2020-03-09 2020-03-10
## 1         31         68         45        140        116         65        376
##   2020-03-11 2020-03-12 2020-03-13 2020-03-14 2020-03-15 2020-03-16 2020-03-17
## 1        322        382        516        548        772       1133       1789
##   2020-03-18 2020-03-19 2020-03-20 2020-03-21 2020-03-22 2020-03-23 2020-03-24
## 1       1362       5964       5526       6327       7676      10567       9893
##   2020-03-25 2020-03-26 2020-03-27 2020-03-28 2020-03-29 2020-03-30 2020-03-31
## 1      12042      18058      17821      19808      19444      20922      26341
##   2020-04-01 2020-04-02 2020-04-03 2020-04-04 2020-04-05 2020-04-06 2020-04-07
## 1      25070      30380      31745      33283      28152      29515      30804
##   2020-04-08 2020-04-09 2020-04-10 2020-04-11 2020-04-12 2020-04-13 2020-04-14
## 1      31533      34126      33755      29861      28917      25306      27051
##   2020-04-15 2020-04-16 2020-04-17 2020-04-18 2020-04-19 2020-04-20 2020-04-21
## 1      28680      31242      32114      32491      26612      25517      27539
##   2020-04-22 2020-04-23 2020-04-24 2020-04-25 2020-04-26
## 1      28486      28819      36188      32796      27631
## 
## $Growth.Rate
##   geo.loc 2020-01-24 2020-01-25 2020-01-26 2020-01-27 2020-01-28 2020-01-29
## 1      US         NA          0         NA          0        NaN        NaN
##   2020-01-30 2020-01-31 2020-02-01 2020-02-02 2020-02-03 2020-02-04 2020-02-05
## 1        NaN         NA        0.5          0         NA          0        NaN
##   2020-02-06 2020-02-07 2020-02-08 2020-02-09 2020-02-10 2020-02-11 2020-02-12
## 1        NaN        NaN        NaN        NaN        NaN         NA          0
##   2020-02-13 2020-02-14 2020-02-15 2020-02-16 2020-02-17 2020-02-18 2020-02-19
## 1         NA          0        NaN        NaN        NaN        NaN        NaN
##   2020-02-20 2020-02-21 2020-02-22 2020-02-23 2020-02-24 2020-02-25 2020-02-26
## 1        NaN         NA          0        NaN         NA          0         NA
##   2020-02-27 2020-02-28 2020-02-29 2020-03-01 2020-03-02 2020-03-03 2020-03-04
## 1  0.1666667          2          4       0.75          4  0.8333333       1.55
##   2020-03-05 2020-03-06 2020-03-07 2020-03-08 2020-03-09 2020-03-10 2020-03-11
## 1   2.193548  0.6617647   3.111111  0.8285714  0.5603448   5.784615   0.856383
##   2020-03-12 2020-03-13 2020-03-14 2020-03-15 2020-03-16 2020-03-17 2020-03-18
## 1   1.186335   1.350785   1.062016   1.408759   1.467617   1.578994  0.7613192
##   2020-03-19 2020-03-20 2020-03-21 2020-03-22 2020-03-23 2020-03-24 2020-03-25
## 1   4.378855  0.9265594   1.144951   1.213213   1.376628  0.9362165   1.217224
##   2020-03-26 2020-03-27 2020-03-28 2020-03-29 2020-03-30 2020-03-31 2020-04-01
## 1   1.499585  0.9868756   1.111498  0.9816236   1.076013    1.25901  0.9517482
##   2020-04-02 2020-04-03 2020-04-04 2020-04-05 2020-04-06 2020-04-07 2020-04-08
## 1   1.211807   1.044931   1.048449  0.8458372   1.048416   1.043673   1.023666
##   2020-04-09 2020-04-10 2020-04-11 2020-04-12 2020-04-13 2020-04-14 2020-04-15
## 1   1.082231  0.9891285  0.8846393  0.9683869  0.8751254   1.068956    1.06022
##   2020-04-16 2020-04-17 2020-04-18 2020-04-19 2020-04-20 2020-04-21 2020-04-22
## 1   1.089331   1.027911   1.011739  0.8190576  0.9588531   1.079241   1.034388
##   2020-04-23 2020-04-24 2020-04-25 2020-04-26 NA
## 1    1.01169   1.255699  0.9062673  0.8425113 NA

3.2 SEIRD Model Simulation in the U.S.

3.2.1 Model Tuning

SEIRD = function(S0, E0, I0, R0, D0, alpha1, alpha2, c1, c2, beta1, gamma1, eta1, n1) {
  Out1 = matrix(0, ncol = 5, nrow = n1) # Output container
  for(i in 1:n1) {
    S0n = S0
    I0n = I0
    R0n = R0
    E0n = E0
    D0n = D0
    ind1 = ifelse(i > c1, 1, 0) # intervetion
    ind2 = ifelse(i > c2, 1, 0) # intervetion
    S0 = max(0, S0n - (alpha1 + alpha2*ind1 + log(i)*0.068*alpha2*ind2)*I0n*S0n)
    E0 = max(0, E0 + (alpha1 + alpha2*ind1 + log(i)*0.068*alpha2*ind2)*I0n*S0n - gamma1*E0n)
    I0 = max(0, I0n + gamma1*E0n - beta1*I0n - eta1*I0n)
    R0 = max(0, R0n + beta1*I0n)
    D0 = max(0, D0n + eta1*I0n)
    Out1[i,1] = S0
    Out1[i,2] = I0
    Out1[i,3] = R0
    Out1[i,4] = E0
    Out1[i,5] = D0
  }
  Out2 = tibble(S = Out1[,1],
                I = Out1[,2],
                R = Out1[,3],
                E = Out1[,4],
                D = Out1[,5])
  return(Out2)
}

# SEIRD Model Parameter Tuning
S0 = 327200000 # Initial susceptible
I0 = 1 # Initial infected
R0 = 0 # Initial recovered
E0 = 10 # Initial exposed
D0 = 0 # Initial death

alpha1 = 1/860000000 # S - > E
alpha2 = -1/1250000000 # Intervention
beta1 = 1/140 # I -> R
gamma1 = 1/7 # E -> I
eta1 = 1/230 # I -> D
n1 = nrow(ts.us.all)

Out1 = SEIRD(S0, E0, I0, R0, D0, alpha1, alpha2, 70, 82, beta1, gamma1, eta1, n1)

plot(1:n1, Out1$I, col = "orange", main = "Infected", xlab="Days", ylab="Cases")
lines(1:n1, ts.us.all$active, col = "orange")

plot(1:n1, Out1$D, col = "red", main = "Deaths", xlab="Days", ylab="Cases")
lines(1:n1, ts.us.all$death, col = "red")

plot(1:n1, Out1$R, col = "green",  main = "Recovered", xlab="Days", ylab="Cases")
lines(1:n1, ts.us.all$recovered, col = "green")

3.2.2 1000 Days Simulation

SEIRD = function(S0, E0, I0, R0, D0, alpha1, alpha2, c1, c2, beta1, gamma1, eta1, n1) {
  Out1 = matrix(0, ncol = 5, nrow = n1) # Output container
  for(i in 1:n1) {
    S0n = S0
    I0n = I0
    R0n = R0
    E0n = E0
    D0n = D0
    ind1 = ifelse(i > c1, 1, 0) # intervetion
    ind2 = ifelse(i > c2, 1, 0) # intervetion
    S0 = max(0, S0n - (alpha1 + alpha2*ind1 + log(i)*0.081*alpha2*ind2)*I0n*S0n)
    E0 = max(0, E0 + (alpha1 + alpha2*ind1 + log(i)*0.081*alpha2*ind2)*I0n*S0n - gamma1*E0n)
    I0 = max(0, I0n + gamma1*E0n - beta1*I0n - eta1*I0n)
    R0 = max(0, R0n + beta1*I0n)
    D0 = max(0, D0n + eta1*I0n)
    Out1[i,1] = S0
    Out1[i,2] = I0
    Out1[i,3] = R0
    Out1[i,4] = E0
    Out1[i,5] = D0
  }
  Out2 = tibble(S = Out1[,1],
                I = Out1[,2],
                R = Out1[,3],
                E = Out1[,4],
                D = Out1[,5])
  return(Out2)
}

# Simulation for 1000 days
n2 = n1+1000
Out1 = SEIRD(S0, E0, I0, R0, D0, alpha1, alpha2, 70, 82, beta1, gamma1, eta1, n2)

plot(1:n2, Out1$I, col = "orange", main = "SEIRD Simulation", xlab="Days", ylab="Cases", ylim = c(0, 1500000), type = "l")
lines(1:n2, Out1$D, col = "red", xlab="Days", ylab="Cases", type = "l")
lines(1:n2, Out1$R, col = "green", xlab="Days", ylab="Cases", type = "l")